Watersheds

2.0 Overview

The following compares two R-based workflows for watershed delineation and hydrological network extraction from digital elevation models (DEMs). We implement parallel processing pipelines using flowdem (Martinsen, 2023) and whitebox (Wu & Brown, 2025) packages to evaluate algorithmic differences in depression handling, flow routing, and stream network generation. Both workflows leverage high-performance geospatial libraries including RichDEM (Barnes, 2016), WhiteboxTools, GDAL 3.11 (contributors, 2020), and the r-spatial ecosystem (Hijmans, 2025; Pebesma & Bivand, 2023a, 2023b).

The comparative approach serves two purposes: (1) validating watershed boundaries through algorithmic consensus, and (2) identifying optimal processing strategies for challenging terrain features such as endorheic basins and artificial depressions. Performance benchmarks and visual comparisons guide selection between rapid exploratory analysis (flowdem) and fine-tuned parametric control (whitebox).

2.1 Declare AOIs

crs_master = sf::st_crs("EPSG:3857")
country = giscoR::gisco_get_countries(
  country = "Malawi", resolution = "3") |>
  sf::st_cast() |> sf::st_transform(crs_master)
lake  = sf::st_read("../assets/inputs/lakes_site.shp") |>
  sf::st_cast() |> sf::st_transform(crs_master)
NA Reading layer `lakes_site' from data source 
NA   `/Users/seamus/repos/map-templates/assets/inputs/lakes_site.shp' 
NA   using driver `ESRI Shapefile'
NA Simple feature collection with 1 feature and 22 fields
NA Geometry type: POLYGON
NA Dimension:     XY
NA Bounding box:  xmin: 35.6 ymin: -15.5 xmax: 35.9 ymax: -15
NA Geodetic CRS:  WGS 84

bbox = lake |>
  sf::st_buffer(dist = 60000) |> 
  sf::st_bbox() |>
  sf::st_as_sfc() |>
  sf::st_sf()

# Interactive map mode: "view"
tmap::tmap_mode("view")
tmap::tm_shape(country) + tmap::tm_borders(lwd=1, col= "green") +
  tmap::tm_shape(bbox) + tmap::tm_borders(lwd=2, col= "orange") +
  tmap::tm_shape(lake) + tmap::tm_borders(lwd=2, col= "blue") +
  tmap::tm_basemap("Esri.WorldImagery")

Figure 2a: Interactive map illustrating layout of area of interest polygons (AOI)

2.2 Download DEM

We acquired elevation data using the elevatr package (Hollister et al., 2023) that accesses via Amazon Web Services Terrain Tiles and the Open Topography Global DEM API collections of global digital elevation models including SRTMGL3, SRTMGL1, AW3D30, and SRTM15Plus. Depending on source, OSM’s Slippy Tiling syntax¹ is used where zoom represents number of vertical tiles counting from top-left corner in a Peudo-Mercator grid (EPSG:3857).

At our study latitude, zoom=11 returns approximately 3-arc-second resolution (~76m at the equator). We then resample to a standardized 100m grid to facilitate downstream metrics computed at landscape scale. The crs_master variable ensures coordinate reference system consistency throughout subsequent operations.

# z = 12: 1-Arc Second Resolution
# z = 11: 3-Arc Second Resolution
# z = 10: 5-Arc Second Resolution
dem = elevatr::get_elev_raster(bbox, z=10, clip="locations") |> terra::rast() 
names(dem) = "elevation" 

dem_100m  = stars::st_warp(
  stars::st_as_stars(dem), 
  cellsize=100, crs=sf::st_crs(crs_master)) |>
  terra::rast()

# Outputs
terra::writeRaster(dem_100m, overwrite = TRUE,
  "../assets/TIF/dem_chilwa_00_raw.tif")

2.3 Hydraulic Conditioning

Both flowdem and whitebox workflows follow the standard hydrological processing sequence: (1) hydraulic conditioning through depression handling, (2) flow direction and accumulation calculation, and (4) stream network extraction. However, they differ fundamentally in algorithmic implementation and parametric control.

flowdem integrates RichDEM operations using priority-flood algorithms (Barnes et al., 2014) with minimal parameters that are optimized for computational efficiency and exploratory analysis. Meanwhile, whitebox provides granular control through extensive tools, enabling parameter tuning for complex terrain such as karst formations or anthropogenic modifications (Hasan et al., 2021).

2.3.1 Depression Breaching

DEMs contain spurious depressions from acquisition artifacts, interpolation errors, and measurement limitations (Lindsay, 2016). These disrupt flow connectivity essential for watershed delineation. Hydraulic conditioning removes artifacts while preserving genuine topographic features like Lake Chilwa’s endorheic basin.

The flowdem workflow implements Barnes et al.’s priority-flood algorithm (Barnes et al., 2014) via the RichDEM library, which processes depressions via optimized queue-based operations. The breach() function carves minimal-depth channels through topographic barriers, improving efficiency with in-memory operations.

The whitebox workflow provides additional parameter controls with which we add a breaching constraint (max_depth=3m) and a flat_increment component to account for flat gradients at lowest elevations. The flat_increment parameter performs similar function as epsilon used in subsequent flowdem::breach operations, only with manual controls over gradient slope(Lindsay, 2016).

# Compare breaching algorithms 
time_wb <- system.time({
  whitebox::wbt_breach_depressions(
    dem     = "dem_chilwa_00_raw.tif",
    output  = "dem_chilwa_11_breached_flat.tif",
    wd      = "../assets/TIF/",
    max_depth = 3,
    flat_increment = 0.001)
})

time_fd <- system.time({
  dem_breach_fd <- flowdem::breach(
    terra::rast("../assets/TIF/dem_chilwa_00_raw.tif"))
  terra::writeRaster(dem_breach_fd, overwrite = TRUE,
    "../assets/TIF/dem_chilwa_01_breached.tif")
})

# Compare runtimes
cat("\n=== Runtime Comparison ===\n")
NA 
NA === Runtime Comparison ===
cat("Whitebox breach:\n")
NA Whitebox breach:
cat(sprintf("  Elapsed time: %.3f seconds\n", time_wb["elapsed"]))
NA   Elapsed time: 1.471 seconds
cat(sprintf("  User time:    %.3f seconds\n", time_wb["user.self"]))
NA   User time:    0.007 seconds
cat(sprintf("  System time:  %.3f seconds\n\n", time_wb["sys.self"]))
NA   System time:  0.005 seconds

cat("flowdem breach:\n")
NA flowdem breach:
cat(sprintf("  Elapsed time: %.3f seconds\n", time_fd["elapsed"]))
NA   Elapsed time: 5.720 seconds
cat(sprintf("  User time:    %.3f seconds\n", time_fd["user.self"]))
NA   User time:    5.383 seconds
cat(sprintf("  System time:  %.3f seconds\n\n", time_fd["sys.self"]))
NA   System time:  0.261 seconds


# Compare spatially
dem_100m   = terra::rast("../assets/TIF/dem_chilwa_00_raw.tif")
dem_breach_fd = terra::rast("../assets/TIF/dem_chilwa_01_breached.tif")
dem_breach_wb = terra::rast("../assets/TIF/dem_chilwa_11_breached_flat.tif")
dem_breach_diff_fd <- dem_100m - dem_breach_fd
dem_breach_diff_wb <- dem_100m - dem_breach_wb
dem_breach_diff_fd[dem_breach_diff_fd == 0] <- NA
dem_breach_diff_wb[dem_breach_diff_wb == 0] <- NA

tmap::tmap_mode("view")
tmap::tm_shape(dem_breach_diff_fd) + tmap::tm_raster(values="Blues", title = "Change (m)",
    breaks = c(-400, -100, -10, 0, 10, 100, 400),
    labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#  tmap::tm_compass(color.dark="gray60,position=c("RIGHT", "top")) +
#  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_title("Breaching Effects (Barnes, 2016)", size=.8) -> tm_breach_fd

tmap::tm_shape(dem_breach_diff_wb) + tmap::tm_raster(values="Blues", title = "Change (m)",
    breaks = c(-400, -100, -10, 0, 10, 100, 400),
    labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +  
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_title("Breaching Effects (Lindsay, 2016a)", size=.8) -> tm_breach_wb

tmap::tmap_arrange(tm_breach_fd, tm_breach_wb, nrow=1)

Figure 2b: Maps comparing effects of breaching algorithms (Barnes, 2016; Lindsay, 2016). Both algorithms target similar depression features but with different modification magnitudes controlled by respective parameters.

2.3.2 Depression Filling

In Lake Chilwa’s endorheic basin, water accumulates via inflow and evaporation without surface outflow. This presents statistical challenges for standard watershed algorithms (O’Callaghan & Mark, 1984), which are typically designed to expect exorheic systems and outlet pour points as inputs (Wu et al., 2015). We adapted our workflow to handle this closed basin by identifying terminal sinks near maximum flow values, enabling internal sub-basin delineation.

• Hydraulic breaching and filling,
• D8 Flow direction and accumulation,
• Depression network characterization,
• Derive terminal drainage points,
• Delineate watershed and streams.

Using the flowdem package, we applied an epsilon-gradient filling via flowdem::fill() to ensure continuous flow through flat areas by adding infinitesimal elevation increments. Meanwhile, the whitebox workflow uses Wang & Liu’s scan-line algorithm (Wang & Liu, 2006) via wbt_fill_depressions_wang_and_liu() to smooth remaining flats in low-relief terrain.

# Compare filling algorithms 
time_wb <- system.time({
  whitebox::wbt_fill_depressions_wang_and_liu(
  dem    = "dem_chilwa_11_breached_flat.tif",
  output = "dem_chilwa_12_filled_wang.tif",
  wd     = "../assets/TIF/")
})

time_fd <- system.time({
  dem_fill_fd <- flowdem::fill(
    terra::rast("../assets/TIF/dem_chilwa_00_raw.tif"), epsilon=T)
  terra::writeRaster(dem_fill_fd, overwrite = T,
    "../assets/TIF/dem_chilwa_02_filled.tif")
})

# Compare runtimes
cat("\n=== Runtime Comparison ===\n")
NA 
NA === Runtime Comparison ===
cat("Whitebox fill:\n")
NA Whitebox fill:
cat(sprintf("  Elapsed time: %.3f seconds\n", time_wb["elapsed"]))
NA   Elapsed time: 1.184 seconds
cat(sprintf("  User time:    %.3f seconds\n", time_wb["user.self"]))
NA   User time:    0.004 seconds
cat(sprintf("  System time:  %.3f seconds\n\n", time_wb["sys.self"]))
NA   System time:  0.004 seconds

cat("flowdem fill:\n")
NA flowdem fill:
cat(sprintf("  Elapsed time: %.3f seconds\n", time_fd["elapsed"]))
NA   Elapsed time: 2.042 seconds
cat(sprintf("  User time:    %.3f seconds\n", time_fd["user.self"]))
NA   User time:    2.017 seconds
cat(sprintf("  System time:  %.3f seconds\n\n", time_fd["sys.self"]))
NA   System time:  0.020 seconds

# Compare spatially
dem_100m   = terra::rast("../assets/TIF/dem_chilwa_00_raw.tif")
dem_fill_fd = terra::rast("../assets/TIF/dem_chilwa_02_filled.tif")
dem_fill_wb = terra::rast("../assets/TIF/dem_chilwa_12_filled_wang.tif")
dem_fill_diff_fd <- dem_100m - dem_fill_fd
dem_fill_diff_wb <- dem_100m - dem_fill_wb
dem_fill_diff_fd[dem_fill_diff_fd == 0] <- NA
dem_fill_diff_wb[dem_fill_diff_wb == 0] <- NA

tmap::tmap_mode("view")
tmap::tm_shape(dem_fill_diff_fd) + tmap::tm_raster(values="RdYlBu", title = "Change (m)",
    breaks = c(-400, -100, -10, 0, 10, 100, 400),
    labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
#  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_title("Filling Effects (Barnes, 2016)", size=.8) -> tm_fill_fd

tmap::tm_shape(dem_fill_diff_wb) + tmap::tm_raster(values="RdYlBu", title = "Change (m)",
    breaks = c(-400, -100, -10, 0, 10, 100, 400),
    labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_title("Filling Effects (Wang & Liu, 2006)", size=.8) -> tm_fill_wb

tmap::tmap_arrange(tm_fill_fd, tm_fill_wb, nrow=1)

Figure 2c: Interactive maps comparing effects of filling algorithms (Lindsay, 2016; Wang & Liu, 2006).

2.4 Flow Accumulation

Flow routing was calculated using the D8 method applied in the flowdem::dirs() and whitebox::wbt_d8_pointer() functions. The D8 single-flow-direction algorithm (O’Callaghan & Mark, 1984) routes flow from each cell to its steepest downslope neighbor among eight adjacent cells. Flow direction pointers are then used to compute flow accumulation values, which represent the number of upstream cells contributing flow to each location. At 100m resolution, each cell represents 0.01 km² (1 hectare) of contributing area, enabling direct conversion between cell counts and drainage area in km².

time_wb <- system.time({
  whitebox::wbt_d8_pointer(
  dem = "dem_chilwa_12_filled_wang.tif",
  output = "dem_chilwa_13_flow_direction_D8.tif",
  wd = "../assets/TIF/")
  whitebox::wbt_d8_flow_accumulation(
  input = "dem_chilwa_13_flow_direction_D8.tif", 
  output= "dem_chilwa_14_flow_accumulation_D8.tif",
  wd    = "../assets/TIF/", pntr=T)
  })

time_fd <- system.time({
  dem_dir_fd <- flowdem::dirs(
    terra::rast("../assets/TIF/dem_chilwa_02_filled.tif"), mode="d8")
  dem_acc_fd <- flowdem::accum(dem_dir_fd, mode="d8")
  terra::writeRaster(dem_dir_fd, overwrite=T,
    "../assets/TIF/dem_chilwa_03_flow_direction.tif")
  terra::writeRaster(dem_acc_fd, overwrite=T, 
  "../assets/TIF/dem_chilwa_04_flow_accumulation.tif")
  })

# Compare runtimes
dem_acc_fd = terra::rast("../assets/TIF/dem_chilwa_04_flow_accumulation.tif")
dem_acc_wb = terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation_D8.tif")

cat("\n=== Runtime Comparison ===\n")
NA 
NA === Runtime Comparison ===
cat("whitebox flow directions:\n")
NA whitebox flow directions:
cat(sprintf("  Elapsed time: %.3f seconds\n", time_wb["elapsed"]))
NA   Elapsed time: 1.248 seconds
cat(sprintf("  User time:    %.3f seconds\n", time_wb["user.self"]))
NA   User time:    0.010 seconds
cat(sprintf("  System time:  %.3f seconds\n\n", time_wb["sys.self"]))
NA   System time:  0.009 seconds

cat("flowdem flow directions:\n")
NA flowdem flow directions:
cat(sprintf("  Elapsed time: %.3f seconds\n", time_fd["elapsed"]))
NA   Elapsed time: 1.792 seconds
cat(sprintf("  User time:    %.3f seconds\n", time_fd["user.self"]))
NA   User time:    1.736 seconds
cat(sprintf("  System time:  %.3f seconds\n\n", time_fd["sys.self"]))
NA   System time:  0.049 seconds

cat("whitebox flow accumulation:\n")
NA whitebox flow accumulation:
cat("Minimum:", min(values(dem_acc_wb), na.rm = TRUE),"cells\n") 
NA Minimum: 1 cells
cat("Median:", median(values(dem_acc_wb), na.rm =TRUE), "cells\n") 
NA Median: 2 cells
cat("Maximum:", max(values(dem_acc_wb), na.rm= TRUE), "cells\n") 
NA Maximum: 1152591 cells
cat("Max contributing area:", round(max(values(dem_acc_wb), na.rm=T) * 0.01, 1), "km²\n")
NA Max contributing area: 11526 km²

cat("flowdem flow accumulation:\n")
NA flowdem flow accumulation:
cat("Minimum:", min(values(dem_acc_fd), na.rm = TRUE),"cells\n") 
NA Minimum: 1 cells
cat("Median:", median(values(dem_acc_fd), na.rm =TRUE), "cells\n") 
NA Median: 1 cells
cat("Maximum:", max(values(dem_acc_fd), na.rm= TRUE), "cells\n") 
NA Maximum: 5017 cells
cat("Max contributing area:", round(max(values(dem_acc_fd), na.rm=T) * 0.01, 1), "km²\n")
NA Max contributing area: 50.2 km²

tmap::tmap_mode("view")
tmap::tm_shape(dem_acc_wb) + tmap::tm_raster(
  values="OrRd", title = "Flow Acc.\n(whitebox)",
  breaks = c(1, 2, 3, 5, 10, 50, 100, 1000, 1200000),
  labels = c("1-2", "2-3", "3-5", "5-10", "10-50", "50-100", "100-1K", ">1K")) +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +  
#  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_acc_wb

tmap::tm_shape(dem_acc_fd) + tmap::tm_raster(
  values="OrRd", title = "Flow Acc. \n(flowdem)",
  breaks = c(1, 2, 3, 5, 10, 50, 100, 1000, 1200000),
  labels = c("1-2", "2-3", "3-5", "5-10", "10-50", "50-100", "100-1K", ">1K")) +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +  
#  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_acc_fd

tmap::tmap_arrange(tm_acc_wb, tm_acc_fd, nrow=1)

Figure 2d: Maps comparing flow accumulation results derived from whitebox (left) and flowdem workflows (right).

2.5 Stochastic Sink Analysis

Given the technical challenges of delineating endorheic watersheds without surface outflows, we apply several additional checks below. Specifically, we attempt to characterize depression structures with complementary topographic metrics to account for DEM uncertainty.

Surface depressions in DEMs contain inherent uncertainty from acquisition artifacts, interpolation errors, and measurement precision limitations. The Stochastic Depression Analysis Tool developed in the whitebox library (Lindsay & Creed, 2005) maps topographic depressions while accounting for uncertainty in depression shape resulting from DEM error (Lindsay & Creed, 2006).

The tool uses a Monte Carlo approach where each grid cell follows a Gaussian error probability distribution function (PDF) with mean of zero and standard deviation equal to the LiDAR RMSE. For our 100m DEM derived from multiple sources, we estimate RMSE conservatively at 3m based on vertical accuracy specifications for SRTM data.

With each iteration, random samples are drawn from the Gaussian error PDF and added to the original DEM. Depressions in the error-added DEM are filled using Wang and Liu’s efficient algorithm (Wang & Liu, 2006), and affected cells are flagged cumulatively. Depressions identified in ≥80% of iterations are classified as “true depressions” rather than artifacts.

We apply 50 iterations following established protocols (Lindsay & Creed, 2005), balancing computational efficiency with detection reliability. This stochastic approach provides confidence estimates for depression features, critical for distinguishing Lake Chilwa’s genuine endorheic basin from spurious artifacts.

# Compute uncertainty using raw and conditioned DEMs
whitebox::wbt_stochastic_depression_analysis(
  dem = "dem_chilwa_00_raw.tif",
  output = "dem_chilwa_15a_depressions_prob.tif",
  wd = "../assets/TIF/",
  rmse = 3.0,
  range = 10.0,
  iterations = 50
)

# Derive binary mask from probability of sinks (~20-80%)
sink_prob = terra::rast("../assets/TIF/dem_chilwa_15a_depressions_prob.tif")
sink_mask_Q1 = sink_prob >= 0.2
sink_mask_Q2 = sink_prob >= 0.4
sink_mask_Q3 = sink_prob >= 0.6
sink_mask_Q4 = sink_prob >= 0.8
sink_mask_Q1[sink_mask_Q1 == 0] <- NA
sink_mask_Q2[sink_mask_Q2 == 0] <- NA
sink_mask_Q3[sink_mask_Q3 == 0] <- NA
sink_mask_Q4[sink_mask_Q4 == 0] <- NA

terra::writeRaster(sink_mask_Q1, overwrite = TRUE,
  "../assets/TIF/dem_chilwa_15a_depressions_Q1.tif")
terra::writeRaster(sink_mask_Q2, overwrite = TRUE,
  "../assets/TIF/dem_chilwa_15a_depressions_Q2.tif")
terra::writeRaster(sink_mask_Q3, overwrite = TRUE,
  "../assets/TIF/dem_chilwa_15a_depressions_Q3.tif")
terra::writeRaster(sink_mask_Q4, overwrite = TRUE,
  "../assets/TIF/dem_chilwa_15a_depressions_Q4.tif")

tmap::tmap_mode("plot")
tmap::tm_shape(sink_mask_Q1) + tmap::tm_raster(
  values="Greens", title = "Predicted Sinks (>0.2)") +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +  
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q1

tmap::tm_shape(sink_mask_Q2) + tmap::tm_raster(
  values="Greens", title = "Predicted Sinks (>0.4)") +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +  
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q2

tmap::tm_shape(sink_mask_Q3) + tmap::tm_raster(
  values="Greens", title = "Predicted Sinks (>0.6)") +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +  
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q3

tmap::tm_shape(sink_mask_Q4) + tmap::tm_raster(
  values="Greens", title = "Predicted Sinks (>0.8)") +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +  
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
  tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q4

tmap::tmap_arrange(
  tm_sink_mask_Q4, 
  tm_sink_mask_Q3, 
  tm_sink_mask_Q2, 
  tm_sink_mask_Q1, 
  nrow=2
  )

Figure 2e: Maps comparing thresholds of sink probabilities derived from Monte Carlo simulation of Gaussian distribution (Lindsay & Creed, 2005)

2.5.1 Sink Depth and Wetness

Building on stochastic depression identification, we quantify two complementary metrics:

• Sink depth - vertical distance from each cell to spillover point
• Topographic wetness index (TWI) - water accumulation potential based on upslope contributing area and local slope

Sink depth uses the raw DEM to preserve actual depression depths below spillover points. TWI calculation uses the conditioned DEM and flow accumulation to ensure topologically consistent routing. Together with stochastic mapping, these metrics characterize spatial extent and hydrological characteristics needed to identify terminal drainage points for watershed delineation.

Tip

When visually checking results in interactive mapping below, try switching a single layer off using top left drop down button in interactive map windows. Toggling between a stack rasters like this enables quick identification of discrepancies derived from candidate algorithms.

# Calulcate depth of sinks
whitebox::wbt_depth_in_sink(
  dem = "dem_chilwa_00_raw.tif",
  output = "dem_chilwa_15_sink_depth.tif", 
  wd = "../assets/TIF/",
  zero_background = F
)

# Calculate topographic wetness index 
whitebox::wbt_wetness_index(
  sca = "dem_chilwa_14_flow_accumulation_D8.tif",
  slope = "dem_chilwa_12_filled_wang.tif",  
  output = "dem_chilwa_16_wetness_index.tif",
  wd = "../assets/TIF/"
)

sink_depth <- terra::rast("../assets/TIF/dem_chilwa_15_sink_depth.tif")
wetness_idx <- terra::rast("../assets/TIF/dem_chilwa_16_wetness_index.tif")
flow_accum <- terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation_D8.tif")

tmap::tmap_mode("view")
tmap::tm_shape(sink_depth) + tmap::tm_raster(values = "steelblue", 
  alpha=1, title = "Depression Sink Depths", style = "quantile", n=7) +
  tmap::tm_shape(wetness_idx) + tmap::tm_raster(values="darkgreen", 
  alpha=1, title="Wetness Index", style = "quantile", n=7) +
  tmap::tm_shape(lake) + tmap::tm_borders(col = "turquoise", lwd = 2) +
  tmap::tm_layout(legend.text.size = 0.8, legend.title.size = 1)

par(mfrow = c(1, 2))
hist(sink_depth[sink_depth <= 25], breaks = 30,
     main = "Depths 0-10m", xlab = "Depth (m)", col = "steelblue")
hist(values(wetness_idx), breaks = 50, xlab = "TWI Value", 
     main = "Topographic Wetness Index", col = "darkgreen")

Figure 2f: Depression sink depths and their surrounding topographic wetness index scores.

2.6 Terminal Pour Points

Endorheic basins require identification of terminal drainage points confirming the ultimate surface flow outlet (Lehner et al., 2022; Zhang et al., 2020) (Li et al., 2025). We locate these adjacent maximum flow accumulation within the lake basin where convergence exceeds 0.5 million cells, confirming watershed-scale drainage convergence.

We test single-point and multi-point configurations to assess delineation sensitivity to outlet placement in flat terrain. Points are digitized interactively using flow accumulation as reference, then saved for workflow reproducibility.

# Inputs
dem_acc = terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation.tif")
streams = dem_acc >= 1000 
streams[streams == 0] <- NA 

# Manually derive terminal pour points
outlets = mapedit::editMap(mapview::mapView(dem_acc)) 
outlets = outlets$all |> # convert to sf
  sf::st_transform(crs_master) |>
  dplyr::select(geometry)
outlets$id <- "terminal_flow"

outlets_add = mapedit::editMap(mapview::mapView(dem_acc)) 
outlets_add = outlets_add$all |> # convert to sf
  sf::st_transform(crs_master) |>
  dplyr::select(geometry)
outlets_add$id <- "terminal_flow"

# Outputs
sf::st_write(outlets, "../assets/SHP/outlets.shp", delete_layer=T, quiet=T)
sf::st_write(outlets_add, "../assets/SHP/outlets_multiple.shp", delete_layer=T)

# Visualize
tmap::tmap_mode("view")
tmap::tm_shape(lake) + tmap::tm_borders(col="royalblue") +
  tmap::tm_shape(outlets) + tmap::tm_symbols(shape="id",fill="orange", size=0.8) +
  tmap::tm_shape(streams) + tmap::tm_raster(values="brewer.reds")

NA Reading layer `outlets_multiple' from data source 
NA   `/Users/seamus/repos/map-templates/assets/SHP/outlets_multiple.shp' 
NA   using driver `ESRI Shapefile'
NA Simple feature collection with 4 features and 1 field
NA Geometry type: POINT
NA Dimension:     XY
NA Bounding box:  xmin: 3970000 ymin: -1740000 xmax: 3990000 ymax: -1720000
NA Projected CRS: WGS 84 / Pseudo-Mercator

Figure 2d: Terminal pour points positioned at maximum flow accumulation within Lake Chilwa basin.

2.7 Watershed Delineation

The flowdem::watershed() function identifies all cells that drain to a specified feature, in this case, the Lake Chilwa boundary polygon. This approach differs from pour-point-based methods by using the entire lake perimeter as the drainage terminus, automatically capturing all flow paths converging on the water body. For endorheic systems where the terminal sink is spatially extensive (Lake Chilwa spans ~60 km N-S), this polygon-based approach may be more robust than single-point methods, as it captures multiple convergent flow paths around the lake’s perimeter.

The resulting watershed polygon represents the complete catchment area draining to Lake Chilwa through surface flow pathways identified by the D8 algorithm. Post-processing converts the raster watershed output to vector format for integration with other spatial datasets and area calculations.

# Inputs
dem_dir = terra::rast("../assets/TIF/dem_chilwa_13_flow_direction_D8.tif")
dem_acc = terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation.tif")

# Tidy
streams = dem_acc >= 1000 
streams[streams == 0] <- NA 
streams = terra::as.lines(streams)

# Compute: Delineate extent of watershed drainage area
watershed = flowdem::watershed(dem_dir, lake) |>
  terra::as.polygons() |>
  sf::st_as_sf()

# Visualize
tmap::tmap_mode("plot")
tmap::tm_shape(watershed) + tmap::tm_borders(col = "blue", lwd = 2) +
  tmap::tm_shape(streams) + tmap::tm_lines(col = "royalblue", lwd = 1) +
  tmap::tm_shape(lake) + tmap::tm_polygons(fill = "royalblue") +
  tmap::tm_layout(legend.text.size = 0.8, legend.title.size = 1) +
  tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
  tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
  tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2)

# Outputs
sf::st_write(watershed, "../assets/SHP/watershed_chilwa_05_flowdem.shp", delete_dsn=T)

Figure 2d: Map illustrating watershed delineation derived in flowdem

2.8 Streams Delineation

# Stream Classification 
whitebox::wbt_extract_streams(
  flow_accum= "dem_chilwa_14_flow_accumulation_D8.tif", 
  output    = "dem_chilwa_18_streams_sensitive.tif",
  wd        = "../assets/TIF/", 
  zero_background = T,
  threshold = 25, 
  )

whitebox::wbt_remove_short_streams(
  d8_pntr   = "dem_chilwa_14_flow_direction_D8.tif",  # flow direction
  streams   = "dem_chilwa_18_streams_sensitive.tif",
  output    = "dem_chilwa_19_streams_d8.tif",
  wd        = "../assets/TIF/",
  min_length= 200,
  )

whitebox::wbt_find_main_stem(
  d8_pntr = "dem_chilwa_14_flow_direction_D8.tif",
  streams = "dem_chilwa_19_streams_d8.tif",
  output = "dem_chilwa_20_streams_trunk.tif",
  wd = "../assets/TIF/"
)

whitebox::wbt_raster_streams_to_vector(
  streams = "dem_chilwa_20_streams_trunk.tif",
  d8_pntr = "dem_chilwa_14_flow_direction_D8.tif",
  output = "../assets/SHP/streams_chilwa.shp",
  wd = "../assets/TIF/"
)

# Inputs
streams_sf = sf::st_read("../assets/SHP/streams_chilwa.shp")
NA Reading layer `streams_chilwa' from data source 
NA   `/Users/seamus/repos/map-templates/assets/SHP/streams_chilwa.shp' 
NA   using driver `ESRI Shapefile'
NA Simple feature collection with 266 features and 2 fields
NA Geometry type: LINESTRING
NA Dimension:     XY
NA Bounding box:  xmin: 3900000 ymin: -1810000 xmax: 4050000 ymax: -1630000
NA CRS:           NA
dem_filled = terra::rast("../assets/TIF/dem_chilwa_12_filled_wang.tif")
sf::st_crs(streams_sf) = 3857

# Visualize
tmap::tm_shape(dem_filled) + tmap::tm_raster(palette = "Greens") +
  tmap::tm_shape(lake) + tmap::tm_borders(col = "turquoise", lwd = 2) +
  tmap::tm_shape(streams_sf) + tmap::tm_lines(col = "steelblue") +
  tmap::tm_layout(legend.text.size = 0.8, legend.title.size = 1)

crs_master = "EPSG:3005"
island_shoreline = sf::st_read("../assets/SHP/island_shoreline.shp", quiet=T)

# z = 12: 1-Arc Second 
# z = 11: 3-Arc Second 
# z = 10: 5-Arc Second 
dem_7arc = elevatr::get_elev_raster(island_shoreline, z=9,clip="locations")|>
  terra::rast()|> 
  terra::crop(terra::vect(island_shoreline)) |>  
  terra::project(crs_master)

#Process DEM by breaching & filling depressions & flats
dem_condt_7arc  = dem_7arc |> 
  flowdem::breach() |>
  flowdem::fill(epsilon=T) 
  #flowdem::fill_basins() # costal basins

# Calculate flow direction & accumulation
dem_dir_7arc = dem_condt_7arc |> flowdem::dirs(mode="d8")
dem_dir_5arc = dem_condt_5arc |> flowdem::dirs(mode="d8")
dem_acc_7arc = dem_dir_7arc |> flowdem::accum(mode="d8")
dem_acc_5arc = dem_dir_5arc |> flowdem::accum(mode="d8")

#tmap::tm_shape(dem) + 
#  tmap::tm_raster(
#    col.scale = tm_scale_continuous(values = "viridis"),
#    col.legend = tm_legend(title = "Elevation (m)", reverse = T)) + 
#  tmap::tm_graticules(lines = T, labels.rot = c(0, 90), lwd = 0.2) +
#  tmap::tm_scalebar(position = c("LEFT", "BOTTOM"), text.size = 0.5) + 
#  tmap::tm_compass(color.dark="gray60",text.color="gray60")+
#  tmap::tm_basemap("Esri.WorldImagery")
writeRaster(dem_condt_7arc, "../assets/TIF/dem_7arc_condt.tif", overwrite=T)
writeRaster(dem_condt_5arc, "../assets/TIF/dem_5arc_condt.tif", overwrite=T)
writeRaster(dem_dir_7arc, "../assets/TIF/dem_7arc_dir.tif", overwrite=T)
writeRaster(dem_dir_5arc, "../assets/TIF/dem_5arc_dir.tif", overwrite=T)
writeRaster(dem_acc_7arc, "../assets/TIF/dem_7arc_acc.tif", overwrite=T)
writeRaster(dem_acc_5arc, "../assets/TIF/dem_5arc_acc.tif", overwrite=T)

2.9 Build 3D map

rivers_strahler = rivers_site |>
  dplyr::mutate(
    width = as.numeric(
      ORD_FLOW
    ),
    width = dplyr::case_when(
      width == 3 ~ 16,
      width == 4 ~ 14,
      width == 5 ~ 12,
      width == 6 ~ 10,
      width == 7 ~ 6,
      TRUE ~ 0
    )
  ) |>
  sf::st_as_sf() |>
  sf::st_transform(crs = "epsg:4326")

h <- nrow(dem_site)
w <- ncol(dem_site)

dem_matrix = rayshader::raster_to_matrix(dem_site)

dem_matrix |>
  rayshader:: height_shade() |>
  rayshader::add_overlay(
    rayshader::generate_line_overlay(
      geometry   = rivers_strahler,
      extent     = dem_site,
      heightmap  = dem_matrix,
      color      = "#387B9C",
      linewidth  = rivers_strahler$width,
      data_column_width = "width"
      ), alphalayer = 1
    ) |>
  rayshader::add_overlay(
    rayshader::generate_line_overlay(
      geometry   = lakes_site,
      extent     = dem_site,
      heightmap  = dem_matrix,
      color      = "#387B9C"
      ), alphalayer = 1
    ) |>
  rayshader::plot_3d(
    dem_matrix,
    zscale       = 14,
    solid        = T,
    shadow       = T,
    shadow_darkness = 2,
    background   = "white",
    windowsize   = c(600, 600),
    zoom         = 0.6,
    phi          = 40,
    theta        = 0 
  )

2.10 Render 3D map

rayshader::render_highquality(
  preview        = T,
  light          = F,
  lightdirection = c(135, 45),
  lightcolor = c("white"),
  lightaltitude = 25,
  ambient_light = 0.1,
  rotate_env     = 0.4,
  intensity_env  = 0.85,
  interactive    = F,
  parallel       = T,
  width          = w,
  height         = h,
  backgroundhigh="#FFFFFF",
  backgroundlow="#FFFFFF"
  )

Figure 4: Three-dimensional map of Lake Chilwa watershed

Environment Setup

pacman::p_load(
  "bslib",
  "cli", "cols4all", "covr", "cowplot",
  "dendextend", "digest", "DiagrammeR", 
  "dtwclust", "downlit",
  "exactextractr", "elevatr",
  "FNN", "future", "flowdem",
  "gdalUtilities", "geojsonsf", "geos", "geodata", 
  "ggplot2", "ggstats","ggspatial", "ggmap", 
  "ggplotify", "ggpubr", "ggrepel", "giscoR",
  "hdf5r", "httr", "httr2", "htmltools",
  "jsonlite",
  "leafem", "leaflet.providers", "libgeos", 
  "luz", "lwgeom", "leaflet", "leafgl",
  "mapedit", "mapview", "maptiles", 
  "methods", "mgcv", "MPI",
  "ncdf4", "nnet",
  "openxlsx",
  "parallel", "plotly", "proj4", "PROJ", "progress", "purrr",
  "randomForest", "rasterVis", "raster", 
  "rayshader", "rayvertex", 
  "RColorBrewer", "rgl", "rmapshaper", "rsconnect", 
  "RStoolbox", "rts", "rgrass",
  "s2", "sf", "scales", "spdep", "stars", 
  "stringr", "supercells",
  "terra", "terrainr", "testthat", "traudem", "taudem", 
  "tidyverse", "tidyterra", "tools",
  "tmap", "tmaptools", "terrainr",
  "whitebox", "xgboost"
  )

# bleeding edge installs
#remotes::install_github("lucarraro/traudem", force=T)
#pak::pkg_install("MPI")
#devtools::install_github("lucarraro/traudem")
#remotes::install_github("opengeos/whiteboxR", build=F)
#remotes::install_github("giswqs/whiteboxR") 
#install.packages("easypackages")

# assign larger memory to tmap rendering
tmap::tmap_options(component.autoscale=F,
  max.raster = c(plot=9500000, view=10000000)
  )

# ----- Whitebox debugging ---- #
# whitebox doesnt like working with files in memory
# so we need to assign path to working directory 
# whitebox::install_whitebox() # install whitebox and GRASS libraries
whitebox::wbt_init() # activate whitebox library in current directory
list.files("/usr/local/taudem/") # Check traudem.pkg installation
# traudem::taudem_sitrep() # Check traudem dependencies
# Ignore above message "Can't find  `MoveOutletsToStrm`"

# ----- TauDEM d8 Check ---- #
# all(sapply(d8, file.exists))
# Expected result: [1] TRUE
# all(file.exists(unlist(d8)))
# Expected result: [1] TRUE
# all(map_lgl(d8, file.exists))
# Expected result: [1] TRUE
# can_register_taudem()
# Expected result: [1] TRUE
# Whitebox loaded tools: `print(wbt_list_tools())`
# Whitebox shortcuts : `print(wbt_help())`
# Whitebox version : `print(wbt_version())`
# Whitebox tool guide: `print(wbt_tool_help("lidar_info"))`\
# Whitebox parameter guide: `print(wbt_tool_parameters("slope"))`
mapviewOptions(fgb = FALSE)
sf::sf_use_s2(use_s2 = FALSE)

Runtime Log

devtools::session_info()
## ─ Session info ───────────────────────────────────────────
##  setting  value
##  version  R version 4.3.0 (2023-04-21)
##  os       macOS 15.7.1
##  system   aarch64, darwin20
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/Vancouver
##  date     2025-11-02
##  pandoc   3.8 @ /opt/local/bin/ (via rmarkdown)
##  quarto   1.7.33 @ /usr/local/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────
##  package           * version    date (UTC) lib source
##  abind             * 1.4-8      2024-09-12 [1] CRAN (R 4.3.3)
##  backports           1.5.0      2024-05-23 [1] CRAN (R 4.3.3)
##  base64enc           0.1-3      2015-07-28 [1] CRAN (R 4.3.3)
##  bit                 4.6.0      2025-03-06 [1] CRAN (R 4.3.3)
##  bit64               4.6.0-1    2025-01-16 [1] CRAN (R 4.3.3)
##  bitops              1.0-9      2024-10-03 [1] CRAN (R 4.3.3)
##  boot                1.3-32     2025-08-29 [1] CRAN (R 4.3.0)
##  brew                1.0-10     2023-12-16 [1] CRAN (R 4.3.3)
##  brio                1.1.5      2024-04-24 [1] CRAN (R 4.3.3)
##  broom               1.0.8      2025-03-28 [1] CRAN (R 4.3.3)
##  bslib             * 0.9.0      2025-01-30 [1] CRAN (R 4.3.3)
##  cachem              1.1.0      2024-05-16 [1] CRAN (R 4.3.3)
##  callr               3.7.6      2024-03-25 [1] CRAN (R 4.3.3)
##  car                 3.1-3      2024-09-27 [1] CRAN (R 4.3.3)
##  carData             3.0-5      2022-01-06 [1] CRAN (R 4.3.3)
##  caret               7.0-1      2024-12-10 [1] CRAN (R 4.3.3)
##  class               7.3-23     2025-01-01 [1] CRAN (R 4.3.3)
##  classInt            0.4-11     2025-01-08 [1] CRAN (R 4.3.3)
##  cli               * 3.6.5      2025-04-23 [1] CRAN (R 4.3.3)
##  clue                0.3-66     2024-11-13 [1] CRAN (R 4.3.3)
##  cluster             2.1.8.1    2025-03-12 [1] CRAN (R 4.3.3)
##  codetools           0.2-20     2024-03-31 [1] CRAN (R 4.3.1)
##  colorspace          2.1-1      2024-07-26 [1] CRAN (R 4.3.3)
##  cols4all          * 0.9        2025-08-28 [1] CRAN (R 4.3.0)
##  coro                1.1.0      2024-11-05 [1] CRAN (R 4.3.3)
##  countrycode         1.6.1      2025-03-31 [1] CRAN (R 4.3.3)
##  covr              * 3.6.4      2023-11-09 [1] CRAN (R 4.3.1)
##  cowplot           * 1.2.0      2025-07-07 [1] CRAN (R 4.3.3)
##  crayon              1.5.3      2024-06-20 [1] CRAN (R 4.3.3)
##  crosstalk           1.2.2      2025-08-26 [1] CRAN (R 4.3.0)
##  curl                7.0.0      2025-08-19 [1] CRAN (R 4.3.0)
##  data.table          1.17.8     2025-07-10 [1] CRAN (R 4.3.3)
##  DBI                 1.2.3      2024-06-02 [1] CRAN (R 4.3.3)
##  deldir              2.0-4      2024-02-28 [1] CRAN (R 4.3.3)
##  dendextend        * 1.19.1     2025-07-15 [1] CRAN (R 4.3.0)
##  devtools            2.4.5      2022-10-11 [1] CRAN (R 4.3.0)
##  DiagrammeR        * 1.0.11     2024-02-02 [1] CRAN (R 4.3.1)
##  dichromat           2.0-0.1    2022-05-02 [1] CRAN (R 4.3.3)
##  digest            * 0.6.37     2024-08-19 [1] CRAN (R 4.3.3)
##  doParallel          1.0.17     2022-02-07 [1] CRAN (R 4.3.3)
##  downlit           * 0.4.4      2024-06-10 [1] CRAN (R 4.3.3)
##  dplyr             * 1.1.4      2023-11-17 [1] CRAN (R 4.3.1)
##  dtw               * 1.23-1     2022-09-19 [1] CRAN (R 4.3.3)
##  dtwclust          * 6.0.0      2024-07-23 [1] CRAN (R 4.3.3)
##  e1071               1.7-16     2024-09-16 [1] CRAN (R 4.3.3)
##  elevatr           * 0.99.0     2023-09-12 [1] CRAN (R 4.3.0)
##  ellipsis            0.3.2      2021-04-29 [1] CRAN (R 4.3.3)
##  evaluate            1.0.5      2025-08-27 [1] CRAN (R 4.3.0)
##  exactextractr     * 0.10.0     2023-09-20 [1] CRAN (R 4.3.1)
##  extrafont           0.19       2023-01-18 [1] CRAN (R 4.3.3)
##  extrafontdb         1.0        2012-06-11 [1] CRAN (R 4.3.3)
##  farver              2.1.2      2024-05-13 [1] CRAN (R 4.3.3)
##  fastmap             1.2.0      2024-05-15 [1] CRAN (R 4.3.3)
##  flexclust           1.5.0      2025-02-28 [1] CRAN (R 4.3.3)
##  flowdem           * 0.2        2025-09-14 [1] Github (KennethTM/flowdem@98cdb20)
##  FNN               * 1.1.4.1    2024-09-22 [1] CRAN (R 4.3.3)
##  forcats           * 1.0.0      2023-01-29 [1] CRAN (R 4.3.0)
##  foreach             1.5.2      2022-02-02 [1] CRAN (R 4.3.3)
##  Formula             1.2-5      2023-02-24 [1] CRAN (R 4.3.3)
##  fs                  1.6.6      2025-04-12 [1] CRAN (R 4.3.3)
##  future            * 1.67.0     2025-07-29 [1] CRAN (R 4.3.0)
##  future.apply        1.20.0     2025-06-06 [1] CRAN (R 4.3.0)
##  gdalUtilities     * 1.2.5      2023-08-10 [1] CRAN (R 4.3.0)
##  generics            0.1.4      2025-05-09 [1] CRAN (R 4.3.3)
##  geodata           * 0.6-2      2024-06-10 [1] CRAN (R 4.3.3)
##  geojsonsf         * 2.0.3      2022-05-30 [1] CRAN (R 4.3.3)
##  geos              * 0.2.4      2023-11-30 [1] CRAN (R 4.3.3)
##  ggmap             * 4.0.1      2025-04-07 [1] CRAN (R 4.3.3)
##  ggplot2           * 3.5.2      2025-04-09 [1] CRAN (R 4.3.3)
##  ggplotify         * 0.1.2      2023-08-09 [1] CRAN (R 4.3.0)
##  ggpubr            * 0.6.1      2025-06-27 [1] CRAN (R 4.3.3)
##  ggrepel           * 0.9.6      2024-09-07 [1] CRAN (R 4.3.3)
##  ggsignif            0.6.4      2022-10-13 [1] CRAN (R 4.3.0)
##  ggspatial         * 1.1.10     2025-08-24 [1] CRAN (R 4.3.0)
##  ggstats           * 0.10.0     2025-07-02 [1] CRAN (R 4.3.3)
##  giscoR            * 0.6.1      2025-08-11 [1] Github (rOpenGov/giscoR@adfed30)
##  globals             0.18.0     2025-05-08 [1] CRAN (R 4.3.0)
##  glue                1.8.0      2024-09-30 [1] CRAN (R 4.3.3)
##  gower               1.0.2      2024-12-17 [1] CRAN (R 4.3.3)
##  gridExtra           2.3        2017-09-09 [1] CRAN (R 4.3.3)
##  gridGraphics        0.5-1      2020-12-13 [1] CRAN (R 4.3.3)
##  gtable              0.3.6      2024-10-25 [1] CRAN (R 4.3.3)
##  hardhat             1.4.2      2025-08-20 [1] CRAN (R 4.3.0)
##  hdf5r             * 1.3.12     2025-01-20 [1] CRAN (R 4.3.3)
##  hexbin              1.28.5     2024-11-13 [1] CRAN (R 4.3.3)
##  hms                 1.1.3      2023-03-21 [1] CRAN (R 4.3.0)
##  htmltools         * 0.5.8.1    2024-04-04 [1] CRAN (R 4.3.3)
##  htmlwidgets         1.6.4      2023-12-06 [1] CRAN (R 4.3.1)
##  httpuv              1.6.16     2025-04-16 [1] CRAN (R 4.3.3)
##  httr              * 1.4.7      2023-08-15 [1] CRAN (R 4.3.0)
##  httr2             * 1.2.1      2025-07-22 [1] CRAN (R 4.3.0)
##  interp              1.1-6      2024-01-26 [1] CRAN (R 4.3.3)
##  ipred               0.9-15     2024-07-18 [1] CRAN (R 4.3.3)
##  iterators           1.0.14     2022-02-05 [1] CRAN (R 4.3.3)
##  jpeg                0.1-11     2025-03-21 [1] CRAN (R 4.3.3)
##  jquerylib           0.1.4      2021-04-26 [1] CRAN (R 4.3.3)
##  jsonlite          * 2.0.0      2025-03-27 [1] CRAN (R 4.3.3)
##  KernSmooth          2.23-26    2025-01-01 [1] CRAN (R 4.3.3)
##  knitr               1.50       2025-03-16 [1] CRAN (R 4.3.3)
##  later               1.4.4      2025-08-27 [1] CRAN (R 4.3.0)
##  lattice           * 0.22-7     2025-04-02 [1] CRAN (R 4.3.3)
##  latticeExtra        0.6-30     2022-07-04 [1] CRAN (R 4.3.3)
##  lava                1.8.1      2025-01-12 [1] CRAN (R 4.3.3)
##  lazyeval            0.2.2      2019-03-15 [1] CRAN (R 4.3.3)
##  leafem            * 0.2.5      2025-08-28 [1] CRAN (R 4.3.0)
##  leafgl            * 0.2.2      2024-11-13 [1] CRAN (R 4.3.3)
##  leaflegend          1.2.1      2024-05-09 [1] CRAN (R 4.3.3)
##  leaflet           * 2.2.2      2024-03-26 [1] CRAN (R 4.3.1)
##  leaflet.providers * 2.0.0      2023-10-17 [1] CRAN (R 4.3.3)
##  leafpop             0.1.0      2021-05-22 [1] CRAN (R 4.3.0)
##  leafsync            0.1.0      2019-03-05 [1] CRAN (R 4.3.0)
##  libgeos           * 3.11.1-3   2025-03-19 [1] CRAN (R 4.3.3)
##  lifecycle           1.0.4      2023-11-07 [1] CRAN (R 4.3.3)
##  listenv             0.9.1      2024-01-29 [1] CRAN (R 4.3.3)
##  logger              0.4.0      2024-10-22 [1] CRAN (R 4.3.3)
##  lubridate         * 1.9.4      2024-12-08 [1] CRAN (R 4.3.3)
##  luz               * 0.5.0      2025-07-29 [1] CRAN (R 4.3.0)
##  lwgeom            * 0.2-14     2024-02-21 [1] CRAN (R 4.3.1)
##  magrittr            2.0.4      2025-09-12 [1] CRAN (R 4.3.0)
##  mapedit           * 0.7.0      2025-04-20 [1] CRAN (R 4.3.3)
##  maptiles          * 0.10.0     2025-05-07 [1] CRAN (R 4.3.3)
##  mapview           * 2.11.2     2023-10-13 [1] CRAN (R 4.3.1)
##  MASS                7.3-60.0.1 2024-01-13 [1] CRAN (R 4.3.1)
##  Matrix              1.6-5      2024-01-11 [1] CRAN (R 4.3.1)
##  memoise             2.0.1      2021-11-26 [1] CRAN (R 4.3.3)
##  mgcv              * 1.9-3      2025-04-04 [1] CRAN (R 4.3.0)
##  microbenchmark      1.5.0      2024-09-04 [1] CRAN (R 4.3.3)
##  mime                0.13       2025-03-17 [1] CRAN (R 4.3.3)
##  miniUI              0.1.2      2025-04-17 [1] CRAN (R 4.3.3)
##  ModelMetrics        1.2.2.2    2020-03-17 [1] CRAN (R 4.3.3)
##  modeltools          0.2-24     2025-05-02 [1] CRAN (R 4.3.3)
##  MPI               * 0.1.0      2022-04-05 [1] CRAN (R 4.3.0)
##  ncdf4             * 1.24       2025-03-25 [1] CRAN (R 4.3.3)
##  nlme              * 3.1-168    2025-03-31 [1] CRAN (R 4.3.3)
##  nnet              * 7.3-20     2025-01-01 [1] CRAN (R 4.3.3)
##  openxlsx          * 4.2.8      2025-01-25 [1] CRAN (R 4.3.3)
##  pacman              0.5.1      2019-03-11 [1] CRAN (R 4.3.3)
##  parallelly          1.45.1     2025-07-24 [1] CRAN (R 4.3.0)
##  pillar              1.11.0     2025-07-04 [1] CRAN (R 4.3.3)
##  pkgbuild            1.4.8      2025-05-26 [1] CRAN (R 4.3.3)
##  pkgconfig           2.0.3      2019-09-22 [1] CRAN (R 4.3.3)
##  pkgload             1.4.0      2024-06-28 [1] CRAN (R 4.3.3)
##  plotly            * 4.11.0     2025-06-19 [1] CRAN (R 4.3.3)
##  plyr                1.8.9      2023-10-02 [1] CRAN (R 4.3.3)
##  png                 0.1-8      2022-11-29 [1] CRAN (R 4.3.3)
##  prettyunits         1.2.0      2023-09-24 [1] CRAN (R 4.3.3)
##  pROC                1.19.0.1   2025-07-31 [1] CRAN (R 4.3.0)
##  processx            3.8.6      2025-02-21 [1] CRAN (R 4.3.3)
##  prodlim             2025.04.28 2025-04-28 [1] CRAN (R 4.3.3)
##  profvis             0.4.0      2024-09-20 [1] CRAN (R 4.3.3)
##  progress          * 1.2.3      2023-12-06 [1] CRAN (R 4.3.1)
##  progressr           0.15.1     2024-11-22 [1] CRAN (R 4.3.3)
##  PROJ              * 0.6.0      2025-04-03 [1] CRAN (R 4.3.3)
##  proj4             * 1.0-15     2025-03-21 [1] CRAN (R 4.3.3)
##  promises            1.3.3      2025-05-29 [1] CRAN (R 4.3.3)
##  proxy             * 0.4-27     2022-06-09 [1] CRAN (R 4.3.3)
##  ps                  1.9.1      2025-04-12 [1] CRAN (R 4.3.3)
##  purrr             * 1.1.0      2025-07-10 [1] CRAN (R 4.3.0)
##  R6                  2.6.1      2025-02-15 [1] CRAN (R 4.3.3)
##  randomForest      * 4.7-1.2    2024-09-22 [1] CRAN (R 4.3.3)
##  rappdirs            0.3.3      2021-01-31 [1] CRAN (R 4.3.3)
##  raster            * 3.6-32     2025-03-28 [1] CRAN (R 4.3.3)
##  rasterVis         * 0.51.6     2023-11-01 [1] CRAN (R 4.3.3)
##  rayshader         * 0.37.3     2024-02-21 [1] CRAN (R 4.3.1)
##  rayvertex         * 0.12.0     2025-02-03 [1] CRAN (R 4.3.3)
##  RColorBrewer      * 1.1-3      2022-04-03 [1] CRAN (R 4.3.3)
##  Rcpp                1.1.0      2025-07-02 [1] CRAN (R 4.3.3)
##  RcppParallel        5.1.11-1   2025-08-27 [1] CRAN (R 4.3.0)
##  RCurl               1.98-1.17  2025-03-22 [1] CRAN (R 4.3.3)
##  readr             * 2.1.5      2024-01-10 [1] CRAN (R 4.3.1)
##  recipes             1.3.1      2025-05-21 [1] CRAN (R 4.3.3)
##  remotes             2.5.0      2024-03-17 [1] CRAN (R 4.3.3)
##  reshape2            1.4.4      2020-04-09 [1] CRAN (R 4.3.0)
##  rex                 1.2.1      2021-11-26 [1] CRAN (R 4.3.3)
##  rgl               * 1.3.24     2025-06-25 [1] CRAN (R 4.3.3)
##  rgrass            * 0.5-2      2025-02-14 [1] CRAN (R 4.3.3)
##  rlang               1.1.6      2025-04-11 [1] CRAN (R 4.3.3)
##  rmapshaper        * 0.5.0      2023-04-11 [1] CRAN (R 4.3.0)
##  rmarkdown           2.29       2024-11-04 [1] CRAN (R 4.3.3)
##  rpart               4.1.24     2025-01-07 [1] CRAN (R 4.3.3)
##  rsconnect         * 1.5.1      2025-08-28 [1] CRAN (R 4.3.0)
##  RSpectra            0.16-2     2024-07-18 [1] CRAN (R 4.3.3)
##  rstatix             0.7.2      2023-02-01 [1] CRAN (R 4.3.0)
##  RStoolbox         * 1.0.2.1    2025-02-03 [1] CRAN (R 4.3.3)
##  rstudioapi          0.17.1     2024-10-22 [1] CRAN (R 4.3.3)
##  rts               * 1.1-14     2023-10-01 [1] CRAN (R 4.3.3)
##  Rttf2pt1            1.3.12     2023-01-22 [1] CRAN (R 4.3.3)
##  s2                * 1.1.9      2025-05-23 [1] CRAN (R 4.3.3)
##  sass                0.4.10     2025-04-11 [1] CRAN (R 4.3.3)
##  satellite           1.0.6      2025-08-21 [1] CRAN (R 4.3.0)
##  scales            * 1.4.0      2025-04-24 [1] CRAN (R 4.3.3)
##  sessioninfo         1.2.3      2025-02-05 [1] CRAN (R 4.3.3)
##  sf                * 1.0-22     2025-08-25 [1] Github (r-spatial/sf@3660edf)
##  shiny               1.11.1     2025-07-03 [1] CRAN (R 4.3.3)
##  shinyjs             2.1.0      2021-12-23 [1] CRAN (R 4.3.0)
##  shinyWidgets        0.9.0      2025-02-21 [1] CRAN (R 4.3.3)
##  sp                * 2.2-0      2025-02-01 [1] CRAN (R 4.3.3)
##  spacesXYZ           1.6-0      2025-06-06 [1] CRAN (R 4.3.3)
##  spData            * 2.3.4      2025-01-08 [1] CRAN (R 4.3.3)
##  spdep             * 1.4-1      2025-08-31 [1] CRAN (R 4.3.0)
##  stars             * 0.6-8      2025-02-01 [1] CRAN (R 4.3.3)
##  stringi             1.8.7      2025-03-27 [1] CRAN (R 4.3.3)
##  stringr           * 1.5.2      2025-09-08 [1] CRAN (R 4.3.0)
##  supercells        * 1.0.0      2024-02-11 [1] CRAN (R 4.3.1)
##  survival            3.8-3      2024-12-17 [1] CRAN (R 4.3.3)
##  svglite             2.2.1      2025-05-12 [1] CRAN (R 4.3.3)
##  systemfonts         1.2.3      2025-04-30 [1] CRAN (R 4.3.3)
##  terra             * 1.8-60     2025-07-21 [1] CRAN (R 4.3.0)
##  terrainr          * 0.7.6      2025-07-25 [1] CRAN (R 4.3.0)
##  testthat          * 3.2.3      2025-01-13 [1] CRAN (R 4.3.3)
##  textshaping         1.0.3      2025-09-02 [1] CRAN (R 4.3.0)
##  tibble            * 3.3.0      2025-06-08 [1] CRAN (R 4.3.3)
##  tidyr             * 1.3.1      2024-01-24 [1] CRAN (R 4.3.1)
##  tidyselect          1.2.1      2024-03-11 [1] CRAN (R 4.3.1)
##  tidyterra         * 0.7.2      2025-04-14 [1] CRAN (R 4.3.3)
##  tidyverse         * 2.0.0      2023-02-22 [1] CRAN (R 4.3.0)
##  timechange          0.3.0      2024-01-18 [1] CRAN (R 4.3.3)
##  timeDate            4041.110   2024-09-22 [1] CRAN (R 4.3.3)
##  tmap              * 4.1        2025-05-12 [1] CRAN (R 4.3.3)
##  tmaptools         * 3.2        2025-01-13 [1] CRAN (R 4.3.3)
##  torch               0.16.0     2025-08-21 [1] CRAN (R 4.3.0)
##  traudem           * 1.0.3      2025-10-28 [1] Github (lucarraro/traudem@a30c9c9)
##  tzdb                0.5.0      2025-03-15 [1] CRAN (R 4.3.3)
##  unifir              0.2.4      2024-02-01 [1] CRAN (R 4.3.3)
##  units               0.8-7      2025-03-11 [1] CRAN (R 4.3.3)
##  urlchecker          1.0.1      2021-11-30 [1] CRAN (R 4.3.3)
##  usethis             3.1.0      2024-11-26 [1] CRAN (R 4.3.3)
##  uuid                1.2-1      2024-07-29 [1] CRAN (R 4.3.3)
##  V8                  6.0.4      2025-06-04 [1] CRAN (R 4.3.3)
##  vctrs               0.6.5      2023-12-01 [1] CRAN (R 4.3.3)
##  viridis             0.6.5      2024-01-29 [1] CRAN (R 4.3.1)
##  viridisLite         0.4.2      2023-05-02 [1] CRAN (R 4.3.3)
##  visNetwork          2.1.2      2022-09-29 [1] CRAN (R 4.3.0)
##  whitebox          * 2.4.3      2025-10-28 [1] Github (opengeos/whiteboxR@294ea4b)
##  withr               3.0.2      2024-10-28 [1] CRAN (R 4.3.3)
##  wk                  0.9.4      2024-10-11 [1] CRAN (R 4.3.3)
##  xfun                0.53       2025-08-19 [1] CRAN (R 4.3.0)
##  xgboost           * 1.7.11.1   2025-05-15 [1] CRAN (R 4.3.3)
##  XML                 3.99-0.18  2025-01-01 [1] CRAN (R 4.3.3)
##  xml2                1.4.0      2025-08-20 [1] CRAN (R 4.3.0)
##  xtable              1.8-4      2019-04-21 [1] CRAN (R 4.3.3)
##  xts               * 0.14.1     2024-10-15 [1] CRAN (R 4.3.3)
##  yaml                2.3.10     2024-07-26 [1] CRAN (R 4.3.3)
##  yulab.utils         0.2.1      2025-08-19 [1] CRAN (R 4.3.0)
##  zeallot             0.2.0      2025-05-27 [1] CRAN (R 4.3.3)
##  zip                 2.3.3      2025-05-13 [1] CRAN (R 4.3.3)
##  zoo               * 1.8-14     2025-04-10 [1] CRAN (R 4.3.3)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
##  * ── Packages attached to the search path.
## 
## ──────────────────────────────────────────────────────────

Barnes, R. (2016). RichDEM: Terrain analysis software. http://github.com/r-barnes/richdem

Barnes, R., Lehman, C., & Mulla, D. (2014). Priority-flood: An optimal depression-filling and watershed-labeling algorithm for digital elevation models. Computers & Geosciences, 62, 117–127.

contributors, G. (2020). GDAL/OGR geospatial data abstraction software library. Open Source Geospatial Foundation. https://gdal.org

Hasan, M. M., Ullah, M. S., Saha, A. K., & Amin, M. M. (2021). Comparing the performances of multiple rainfall-runoff models of a karst watershed. Asian-Australasian Journal of Bioscience and Biotechnology, 6(1), 26–39.

Hijmans, R. J. (2025). Terra: Spatial data analysis. https://CRAN.R-project.org/package=terra

Hollister, J., Shah, T., Nowosad, J., Robitaille, A. L., Beck, M. W., & Johnson, M. (2023). Elevatr: Access elevation data from various APIs. https://doi.org/10.5281/zenodo.8335450

Lehner, B., Messager, M. L., Korver, M. C., & Linke, S. (2022). Global hydro-environmental lake characteristics at high spatial resolution. Scientific Data, 9(1), 351.

Li, Y., Huang, C., Wang, N., Chen, Y., Zhang, S., & Duan, H. (2025). Enhanced understanding of hydrological connectivity among lakes on the qiangtang plateau. International Journal of Digital Earth, 18(1), 2450845.

Lindsay, J. B. (2016). Efficient hybrid breaching-filling sink removal methods for flow path enforcement in digital elevation models. Hydrological Processes, 30(6), 846–857.

Lindsay, J. B., & Creed, I. F. (2005). Removal of artifact depressions from digital elevation models: Towards a minimum impact approach. Hydrological Processes: An International Journal, 19(16), 3113–3126.

Lindsay, J. B., & Creed, I. F. (2006). Distinguishing actual and artefact depressions in digital elevation data. Computers & Geosciences, 32(8), 1192–1204.

Martinsen, K. T. (2023). Flowdem: Flow routing on digital elevation models.

O’Callaghan, J. F., & Mark, D. M. (1984). The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing, 28(3), 323–344.

Pebesma, E., & Bivand, R. (2023a). Spatial data science: With applications in r. https://doi.org/10.1201/9780429459016

Pebesma, E., & Bivand, R. (2023b). Spatial data science: With applications in r. 352. https://doi.org/10.1201/9780429459016

Wang, L., & Liu, H. (2006). An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. International Journal of Geographical Information Science, 20(2), 193–213.

Wu, Q., & Brown, A. (2025). ’Whitebox’: ’WhiteboxTools’ r frontend. https://CRAN.R-project.org/package=whitebox

Wu, Q., Liu, H., Wang, S., Yu, B., Beck, R., & Hinkel, K. (2015). A localized contour tree method for deriving geometric and topological properties of complex surface depressions based on high-resolution topographical data. International Journal of Geographical Information Science, 29(12), 2041–2060.

Zhang, J., Condon, L. E., Tran, H., & Maxwell, R. M. (2020). A national topographic dataset for hydrological modeling over contiguous united states. Earth System Science Data Discussions, 2020, 1–26.

---

1. https://wiki.openstreetmap.org/wiki/Slippy_map